Eigenvalue distributions from a star product 
approach 



J N Kriel 1 and F G Scholtz 1 2 

1 Institute of Theoretical Physics, University of Stellenbosch, Stellenbosch 7600, 
South Africa 

2 National Institute for Theoretical Physics (NITheP), 7600 Stellenbosch, South 
Africa 

Abstract. We use the well-known isomorphism between operator algebras and 
function spaces equipped with a star product to study the asymptotic properties of 
certain matrix sequences in which the matrix dimension D tends to infinity. Our 
approach is based on the su(2) coherent states which allow for a systematic 1/D 
expansion of the star product. This produces a trace formula for functions of the 
matrix sequence elements in the large-!? limit which includes higher order (finitc- 
D) corrections. From this a variety of analytic results pertaining to the asymptotic 
properties of the density of states, eigenstates and expectation values associated with 
the matrix sequence follows. It is shown how new and existing results in the settings of 
collective spin systems and orthogonal polynomial sequences can be readily obtained 
as special cases. In particular, this approach allows for the calculation of higher order 
corrections to the zero distributions of a large class of orthogonal polynomials. 



1. Introduction 

The eigenvalue analysis of operators that arise in physical problems is often simplified 
by exploiting special structures present in the operator's matrix representation. One 
such structure encountered in some matrix sequences {Md}d =1 , where is D x D 
dimensional, concerns the limiting behaviour of the matrix elements [Mr)} n ,n+m as n 
and D tend to infinity in a fixed ratio. In particular, there may exist a set of functions 
{fm(x)}m=o such that [Mr)\n,n+m —> fm(x) when n, D — t- oo and n/D — > x. Intuitively, 
this suggests that the matrix elements on the m'th diagonal represent an increasingly 
fine sampling of f m (x) at the points x n = n/D. If each f m (x) is continuous the matrix 
elements therefore exhibit a "smooth" structure in that [Mn] n ,n+m is a slowly varying 
function of n when D is large. Since the set {f m {x)} characterises the matrix sequence 
as D — )• oo one expects that the asymptotic behaviour of the eigenvalue distribution, 
eigenstates and expectation values can be studied directly in terms of these functions. 
The development of such a formalism is the main goal of this paper. 
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Despite their usefulness these ideas appear to have found only limited application 
in physical problems. Hollenberg and Witte [1] made use of this structure in the Hamil- 
tonian matrices of extensive many-body systems to derive an exact analytic expression 
for the ground state energy density by minimizing a particular linear combination of 
the / m (a;)'s. Deift and Mclaughlin [2] (see also [3]) studied a continuum limit of the 
Toda lattice system of which the equations of motion amount to the isospectral flow 
of a particular tridiagonal matrix with a dimension given by the system size. Under 
the conditions outlined above the continuum limit of this matrix equation is a partial 
differential equation for the / m (x)'s. In fact, from the matrix viewpoint the Toda lattice 
dynamics represent a special case of the flow equation renormalization scheme proposed 
by Wegner j3] which has been applied to a wide range of quantum mechanical many-body 
problems [5] . Here too the presence of a smooth structure in the matrix elements allows 
the matrix (or operator) flow to be recast as a partial differential equation in the ther- 
modynamic limit [6j [7J . An important feature of these studies is their non-perturbative 
nature, i.e. no assumptions regarding the magnitude of coupling constants are made 
and expansions are instead controlled by the small parameter 1/D. These techniques 
are therefore particularly well suited to the study of quantum critical behaviour in the 
thermodynamic limit [7]. 

Matrix sequences of this type also arise naturally in the theory of orthogonal poly- 
nomials. Kuijlaars and van Assche [SI [9] made use of this structure in the Jacobi matrix 
to derive an expression for the asymptotic zero distribution of certain orthogonal poly- 
nomial sequences. See also Bourget [10] for a recent application to the theory of Jacobi 
matrices of which the elements satisfy a similar small deviation condition. 

In this paper we present a treatment of matrix sequences which exhibit this smooth 
structure using a coherent state star-product formalism. The advantage of this approach 
is that it produces, in a simple and direct manner, analytic results for the asymptotic 
eigenvalue distribution (i.e. the density of states) as well as for eigenstates and expecta- 
tion values. Indeed, our approach is based on a systematic 1/D-expansion which allows 
for the inclusion of finite-/} corrections. This is in contrast to most existing results which 
only consider the D — > oo limit. In particular we are able to extend the results of [8] 
and derive corrections to the asymptotic zero distributions of a large class of orthogonal 
polynomial sequences. In quantum mechanical applications this allows for the inclusion 
of finite-size corrections in the density of states and expectation values. 

The paper is organized as follows. In section [2] we summarise the basic properties 
of spin coherent states and introduce the notions of symbols and the star-product. The 
precise class of real tridiagonal matrix sequences under consideration is defined in section 
|3] and asymptotic expressions for their symbols are derived. Section @] is dedicated to 
the derivation of a trace formula which is used in section [5] to derive our main result: 
an expression for the asymptotic eigenvalue distribution which includes higher order 
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corrections. Section E] deals with expectation values and the structure of eigenstates. 
Results obtained for tridiagonal sequences can be generalized to sequences of banded 
Hamiltonian matrices satisfying certain boundary conditions. We show how this is done 
in section [7J Sections [8] and [9] present applications of our results to collective quantum 
spin systems and orthogonal polynomials. 



2. Background 

2.1. Spin coherent states 

Let Hj be a D = 2j '• + 1 dimensional Hilbert space carrying the j-irrep of su(2) with 
generators { J±, J z }. The spin coherent states [H] [12] are defined as 

\z) = (1 + zz)-ieMzJ+\\j, -3) z e C. (1) 

Here z is the conjugate of z and {\j,m)y m= _j is the standard basis of J z eigenstates. 
The inner product of any two such states is 

\( z \ v )\* = (1 + Z vfi(l + vz) 2j (l + zz)- 2] {\ + vv)~ 2j . (2) 

The set of coherent states forms an over-complete basis for %j and provides a resolution 
of the identity in the form / = J dz dz fi(z, z)\z)(z\ with /j,(z, z) = (2j + l)7r _1 (l + zz)~ 2 . 

Let z = rexp(i8) and set x = r 2 /(l + r 2 ) G [0, 1]. The expansion coefficients of \z) 
in the J 2 -basis are 

c n = 0', —j + n\z) = x n ' 2 {\ - x) i2 ^/ 2 e me (^j ^ n = 0, 1, 2, . . . , 2j. (3) 

Note that P(n; 2j, x) = \c n \ 2 amounts to a binomial distribution with success probability 
x and 2j trails. It is known that if x e (0, 1) is kept fixed while j tends to infinity then 
P(n; 2j, x) will approach a normal distribution with mean 2jx and variance 2jx(l — x). 

2.2. Symbols and the star product 

With each operator H acting on Tij we associate a function H(z,z) = (z\H\z) known 
as the symbol of H. The coherent state resolution of the identity allows the trace of H 
to be expressed in terms of its symbol as 

tr(#) = [ dzdz/i(z,z)H(z,z) = [ dx I dOH{x,9). (4) 

J 2tt J J 

In what follows we use the (z, z) and (x, 8) parametrizations interchangeably and will 
often suppress the functional dependence of the symbols. 

The operator product is realised on the symbol level through the star product 

(z\AB\z) = (z\A\z) * (z\B\z) = A* B (5) 
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which can be expressed in terms of differential operators acting on the symbols as 

* = J dv dv exp[t^ z ] n(v + z, v + z)\(z\v + z)\ 2 exp[v~$z\. (6) 

This exact representation can be brought into a more practical form by noting that 
\(z\v + z)\ 2 is maximal at v — and that expanding its logarithm around this point 
produces, to lowest order, the quadratic expression 

-2j 



log \(z\v + z) 



:VV + 



(7) 



(1 + zzf 

This suggests that a saddle-point approximation of the integral in (jBJ) should allow the 
star product to be expressed as a power series in 1/j. The result can be written as 

1 + zzT+ m tSSJ™ (8) 

k=l J n,m=l 

where is a k x k Hermitian matrix which generally contains z and z. The first three 
of these are A (1) = 1/2, 

~ 12zz(3zz+l) Qz(Qzz + l) 6z 2 

and A (3) 



1+ Ei E 

i — 1 J ^ ™ — 1 



A (2) 



4zz 2z 
2z 1 



1 

48 



Qz(Qzz + 1) 
Qz 2 



3Qzz + 3 
6z 



6z 
1 



(9) 



To linear order in 1/j the star product is therefore given by 



<_ (1 + ^)2^ 



(10) 



2j 

^. 3. Smooth symbols and fluctuations 

Consider a sequence of operators {A/I^ci where Aj acts on the D = 2j '• + 1 dimensional 
space "Hj. We assume that the corresponding sequence of symbols {Aj} as well as their 
derivatives to z and z scale intensively with j, i.e. as O(j ). This smoothness condition 
allows the symbol (z|yl||2:) = Aj * Aj to be approximated to linear order in 1/j as 

! l + zz) 2 



Aj * L Aj 



A 2 + 



2j 



{d z Aj) {d- z Aj 



(11) 



As shown next this result has important consequences for the scaling behaviour of fluc- 
tuations. Since the symbol Aj is an expectation value with respect to \z) we write 
5Aj = Aj - Aj, in which case (z\(dAj) 2 \z} = (1 + zzf (d z Aj) (d- z Aj) /(2j) + 0(j- 2 ). 
The fluctuations in Aj with respect to the coherent states are therefore suppressed by 
a factor of 1/j. This result may well have been anticipated due to the semi-classical 
nature of these states, but here we wish to emphasise the importance of the smoothness 
condition. It follows from induction that higher-order fluctuations are even further sup- 
pressed, with (z\(6Aj) n \z) generally scaling like 0(j~^ n ^ 2 ^). 
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We arrive at the useful result that for a smooth function f(x) the symbol of f{A) 
is given, up to linear order in 1/j, by 

M 4- zz) 2 

(z\f(Aj)\z) = f{Aj) + ±Lp-f{Aj) {d z Aj) (d- z A 3 ) + O^). (12) 
This simple expression forms the cornerstone of the asymptotic expansions that follow. 

3. Matrices with smoothly varying elements 

Next we define the class of matrix sequences to be considered in the sequel. For simpli- 
city we will focus mainly on the tridiagonal case, although all the results can be adapted 
to apply to band diagonal Hermitian matrices of which the off-diagonal elements sat- 
isfy a particular constraint. This straightforward generalization is presented in section [7J 

Let {Hj} denote a sequence of real, symmetric tridiagonal matrices with Hj acting 
on the D = 2j + l dimensional space Hj. The basis in which Hj is given is identified with 
the basis { \j, m) } J m= _j of J z eigenstates. The diagonal matrix elements of Hj are denoted 
by a n j = [Hj] n+ x )n+ i with n = 0, . . . , 2j while the off-diagonal ones are b n j = [Hj] n+1>n 
with n — 1, . . . , 2j. We require that these matrix elements satisfy a certain smoothness 
condition in that they are slowly varying functions of n. To be precise, we assume that 
there exists smooth functions A , A\, B , B\ : [0, 1] — > R such that, up to linear order in 

the matrix elements are given by 

where B (x) > for x e (0, 1). Note that Hj scales intensively with j. 

Our first task is to find a large- j expansion for the symbol Hj of Hj. Some algebra 
reveals that Hj is given exactly by 

1/2 



{z\H j \z) = H j (x,9) = (a n ,j) x + 2 



X 


1/2 / 


~2j - n 


1 — X 




71 + 1 



b n+1J cos(^) (14) 



where {f(n)) x denotes the expectation value of f(n) over n with respect to the binomial 
distribution P(n;2j, x) = \c n \ 2 . Now consider a fixed x G (0,1). At large j the 
distribution P(n; 2j, x) is sharply peaked around the mean (n) = 2jx with fluctuations 
in n/2j of order (9(j -1 / 2 ). This suggests that the expectation values in fTl4j) can be 
approximated using an expansion in orders of 5n = n — (n). Up to linear order in 1/j 
and for x G (0, 1) this produces 

Hj(x,6) =HW{x,0) +H^{x,6)/j + 0{j- 2 ) (15) 

where 

# (0) (x, 9) = A (x) + 2^0(0;) cos(^) 

H®(x, 9) = A x {x) + 2B l {x) cos(9) + %H®{x, 9) - B ^°< 6 \ (16) 
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ff/x,0) 




(a) (b) 

Figure 1. Results for the Toeplitz sequence with a = and b = 1 at finite j. (a) 
The symbol Hj(x, 9) as a function of x at 9 = 0. (b) The logarithm of (z\(6Hj) 2 \z) for 

Z = y/x/(l- X). 



Note that the functions H^(x, 9) and H^\x, 9) can also be defined for x £ [0, 1] through 
the limits 

H(°\x,6) = lim HJx,G) and ffWfo 9) = lim j(i2"ifo 0) - H^(x, 9)). (17) 

j—*oo j—¥00 

A natural question is whether the results of (fTEj) . which were valid for a; £ (0, 1), also 
hold at x = 0, 1. We will show that this is not generally the case. That the end-points 
require special attention is a consequence of the fundamentally different nature of the 
coherent state \z) for x = 0, 1 compared to x £ (0, 1). When x £ (0, 1) the state \z) is 
a linear combination of all 2j + 1 basis states and the symbol Hj = (z\Hj\z) therefore 
depends on all the matrix elements of Hj. In contrast, at x — 0, 1 the state \z) is a 
single basis state, either \j, —j) or \j, +j), and (z\Hj\z) is simply the top left or bottom 
right diagonal matrix element. It follows that 

H(°\0,9) = A (0) H^(1,9) = A (1) 

H^ 1 \0 1 9) = A l (0) H ( < 1 \1,9) = A 1 (1) (18) 

which is to be compared with the results for x £ (0,1) in ( fl6|) . The continuity of 
(x, 9) at x = 0, 1 clearly depends on the behaviour of the functions B (x) and 
Bi(x) at the edges. In particular, continuity to lowest order requires that B (x) vanishes 
at x = 0, 1. We say that a matrix sequence is closed if -Bq(O) = -Bq(I) = 0. This implies 
that the off-diagonal matrix elements approach zero at the top left and bottom right 
corners of the matrix as j — > oo. Continuity in the higher order terms will clearly place 
further restrictions on the behaviour of B (x) and Bi(x) at x = 0, 1. For our purposes 
the notion of closure at lowest order in 1/j is sufficient. 

To illustrate how the continuity of the symbols depend on the behaviour of B (x) at 
x = 0, 1 we consider as an example a sequence of Toeplitz matrices. Since these issues 
are local in nature what is observed here is quite generic and common to all non-closed 
matrix sequences. For tridiagonal Toeplitz matrices with constant entries and 



Eigenvalue distributions from a star product approach 



7 



b n j = b > we have H^(0,9) = H^(l,9) = a while H^(x,9) = a + 2bcos(9) when 
x G (0, 1). The sequence of continuous symbols Hj therefore converge to a function 
which is continuous everywhere except at x = 0, 1. This is reflected by the numerical 
results in figure [1] (a). Furthermore, in a neighbourhood of x = 0, 1 the symbols Hj 
clearly vary very rapidly and have derivatives which scale like j, despite the symbols 
themselves being intensive. From the discussion in section [2731 this is expected to impact 
on the scaling behaviour of the fluctuations close to the edges. This is indeed the case, 
as can be seen in figured] (b). The expansion in ( fl2i) for (z\f(Hj) \z) therefore fails at 
x = 0,1. 



4. Asymptotic expansion of traces 

Let {Hj} be a matrix sequence as defined section[3]and f(x) a smooth bounded function. 
We are interested in the large- j behaviour of the scaled trace 



2j 



TT tr (/(4)) = ^: [to 

+ 1 2n In 



2- 



d9{z\f(Hj)\z) 



which is captured by the expansion 

7- = r (0) +r (1) /j + o(r 2 ) 

where 



T (0) = lim 75 



and 



j-s.oo 



r (1) = iim j(Tj - r (0) ). 



(19) 

(20) 
(21) 

(22) 



For x G (0, 1) the expansion in ffl~2l produces 

(z\f(Hj)\z) = f(H^) + ± \f{H®)HV> + \f"{H^){l + zzf (d.HM) (d- z H^) 

with H(°> and H^ given by ffl6|) . We will show that for closed sequences this result 
suffices to calculate both and T^. However, for sequences which are not closed 
the uncontrolled fluctuations at x — 0, 1 result in non-trivial edge corrections to 
which are not captured by the expansion above. 



4-1. The lowest order contribution 

To determine we first split the integral in (fl9l) into three parts to isolate the 
problematic edges at x = 0,1 from the interior where (I22p holds: 



jio) 



lim 

Jo 



"dx r ^-(z\f(Hj)\z) 



o 



2tt 



lim lim 

e— >-0+ j—too 



dx(. . .) + 



o 



l-e 



dx(...) 



dx(...) 



l-e 



lim lim / dx (...) + lim lim 

e->0+ j— >oo Jq e— 5>0+ J->oo 

S V ' V 



l-e pi 

dx (...)+ lim lim / dx (. . .) 

e^0+ j^oo J x _ t 



(23) 
(24) 
(25) 



r (0 



V 



(0) 



T-(0) 
' R 
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Here , and denote the contributions to coming from the interior and 
left and right edges of the x G [0, 1] interval respectively. Since / is bounded the two 
edge terms are zero. From the expressions for #(°> and in (PJ we also see that 
the 0(j~ l ) terms in ( )22l) can only diverge (if at all) at x — 0, 1 and may therefore be 
neglected when taking the j — > oo limit of J e dx (. . .) for any finite e. It follows that 

pi i-2ir j/i pi p2ir -in 

T{0) = l dX Jo h f{Hm{x ' 6)) = J, dX J ^f(Mx)+2B (x)cos(6)). (26) 

To lowest order in 1/j the edge effects are therefore irrelevant and fl26|) above holds 
regardless of whether or not the matrix sequence is closed. 



4-2. The linear order correction T*- 1 * 1 

From (I2TI) and (|26|) the linear order correction to Tj can be written as 



I v = hm ax — j 
,oc Jo Jo 



2tt' 



(zlfiH^-fiH 



(o)^ 



Here we again isolate the edge corrections by writing = 7^ + 7} 
done for T^ -* in ( )25l) . We consider these three contributions individually. 



(27) 



-(i) 



Tft ' as was 



4-2.1. The interior term Tj Since the integral for Tj L> excludes the edges it follows 
from ([22]) that 



-(i) 



V 



(1) 



lim 

e->0+ 



1-6 ,2* dQ 

ax / — 

A 27T 



f\H^)H^ + ~r(# (o) ) 



(<%g (0) ) 2 

8x(l — 



+ 



where the j-dependence has dropped out. At this stage the integrand contains terms 
(one inside H^) that diverge at the edges and apparently prevent the e — > + limit 
from being taken. However, after applying integration by parts to the ^-integral of the 
second term these divergences are found to cancel. Taking the e — > + limit then yields 



(1) 



1 

dx 
o Jo 



2tt 



d9_ 

2^ 



f( #(0)) 



-{2x - l){d x H^) + A x (x) + 25 1 (x) cos(fl) 



4 



(29) 



4-2.2. The edge terms andT^p The smoothness condition contained in f|T3|) allows 
the edge corrections to be calculated using a simple combinatoric argument which is 
presented in the appendix. It is found that these corrections depend only on the values 
of A (x) and B (x) at x = 0,1. In particular, B (0) = (-Bo(i) = 0) implies that 
Tjp = (7^=0) and the edge corrections therefore vanish if the sequence is closed. 
We define next 



a(x) = An(x) - 2B n (x) 



and (5(x) = A Q (x) + 2B (x) 



and 



r K/3](A) 



ti 1 arcsin 



a+/3-2\ 
a-/3 



A G [a, (3} 
otherwise 



(30) 



(31) 
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If £(,(0) > the correction can be expressed compactly as 

/+00 ^ 
dxf{x)—T [a{0) , m] (x) (32) 

where it should be understood that the discontinuities of T[ a (o),p(o)] (x) at x — a(0),/3(0) 
produce (^-functions in the integrand. If B (l) > the correction 7^ follows from f[3"2l 
by replacing [«(0), (3(0)} by [a(l), 0(1)]. 

This concludes the derivation of the trace formula corresponding to the expansion 
of Tj in (1201) . In the next section we show how detailed information regarding the 
asymptotic eigenvalue distribution of the matrix sequence follows from this result. 

5. The density of states 

Let {A(n, j)}^ 1 denote the eigenvalues of Hj. The scaled density of states associated 
with Hj is defined as 



p J (X) = ^-jJ2S(X-X(n,j)). (33) 

n=l 

Our goal is to derive a linear order large- j expansion of the form 

p 1 (\) = p^(X)+p^(X)/j + 0(r 2 ) (34) 
where p^(A) and p^'(X) are defined by requiring that 



lim / dXp,(X)g(X) = [ dXp(°\X)g(X) 



(35) 



lim / dXj[pj(X) - p <f»(\)]g(\) = / dXpV(X)g(X) (36) 



holds for continuous g{X). It is useful to introduce the distribution functions associated 
with these asymptotic densities as 

X>(°)(A) = [ X dX'p {0) (X') and T? (1) (A) = t dA'p (1) (A'). (37) 



The traces considered in section H] may now be expressed as 

7- = ^i^trU^-)) = J dX Pj (X)f(X) and so T^ l) = j dXp^\X)f{X). (38) 

Knowledge of Z>( 0,1 )(A) is therefore sufficient to calculate 7 - *- ' 1 ) for any choice of f(x). 
To determine P( 0,1 )(A) we first use this relation in reverse and express the density of 
states in terms of a particular trace, namely 

7}(A, 7 ) = 2j^T tr K( A -^)] (39) 

where 7 (x) = (1 + tanh(x/7))/2. Note that 6 7 (x) is smooth for 7 > but converges 
to the Heaviside step function Q(x) as 7 — > + . From fl37|) and fl38l it follows that 

P (0 ' 1} (A) = lim /"dA'p (0 ' 1) (A / )e 7 (A- A') = lim r (0,1) (A,7) (40) 

7-i>0+ J 7^0+ 
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and so from (|2"6"]) we find 

P(°)(A)= fdx r^-Q(\- H (°\x,9)). 
Jo Jo 271 

Performing the ^-integral above yields the final form of T)^(X) as 

£> (0) (A) = / dxU[ a ( x ),i3( x )](\) 



(41) 



(42) 



with a(x) and (3(x) as defined in (|3~0]) and 




7T 1 arcsin 



Q+/3-2A 
a-/3 



A < a 
X> (3 

otherwise 



(43) 



From this it is evident that 

A_ = min{a(x) : x G [0, 1]} and A + = m&x{(3(x) : x G [0, 1]} (44) 

defines the range of eigenvalues in the j — > oo limit. In particular, T>(°>(\) = for 
A < A_ while V^(X) = 1 at A > A+. Expression flU was obtained by Kuijlaars and 
van Assche in [8]. In fact, this result is far more robust than the present derivation sug- 
gests; in [9] it was shown that (142]) also holds for a wide class of discontinuous a(x) and 
(3(x). The expression for A_ as a ground state energy density was previously obtained 



in 



EH- 



To treat the linear order correction T)^(X) we first decompose it as 



P (1) (A) 



lim 

7^0+ 



T L (1) (A,7) + 77U7) + rr(A, 7 ) = V^(X) + V?>(X) + V£(X).(45) 



The interior contribution then follows from (|29j) as 



V^iX) = lim 7} (1) 

7->0+ 



— dx 



2- 



d9_ 
2^ 



5(A-if (0) (a;,^))F(x,cos(^)) 



dx n [a{x)Mx)] (A) F (x, (A - A (x))/(2B Q (x))) 

with F(x, cos{6)) = (2x - l){d x H^)/A + A x (x) + 2£i(x) cos(0). Here 

Ti- l [{X-a){f3-X)Y 1/2 X G (a, f3) 



^[a„fl](A) 



otherwise 



(46) 
(47) 

(48) 



when a < /3 while ft[ a , a ](A) = 5(X - a). Note that V^\X) is zero outside the interval 
[A_,A+]. It only remains to determine the two edge corrections to T)^(X) which are 
given by the 7 — > + limits of 7^ (A, 7) and 7^ (A, 7) for the trace in (J3HD • It follows 
from ([32]) that 



2?«(A) 



47T 



arcsm 



A - A (0) 
2£? (0) 



and Pi (A) = — arcsin 

4:71 



A -4,(1) 
2S (1) 



(49) 
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for A G [a(0),/3(0)] and A G [a(l),/3(l)] respectively. Outside these ranges both correc- 
tions are zero. Also, if B (0) = (A)(1) = 0) then £>[ 1} (A) (V%\\)) is identically zero. 



To summarize, equations f|42|) . fj45|) . f|47j) and (149]) together yield expressions for 
X>( 0,1) (A). Taking derivat ives to A produces /^'^(A) and by using (1381) the corrections 
7~(o,i) f or an arbitrary trace can be calculated. In this regard, note that the linear order 
correction D^^A) is generally not continuous at A = A± and that this gives rise to 
^-functions in the density of states. 



6. Further results for special cases 

In this section we consider matrix sequences for which certain quantities associated 
with the eigenstates are well-behaved functions of the corresponding eigenvalue or state 
label. The conditions under which these requirements are met will be investigated in 
the context of a specific example in section 16.31 



6.1. The index function 

If the eigenvalues of Hj satisfy X(n,j) < X(n + 1, j) we define the index function as 

l(X(n,j),j) = ^ n , (50) 

which is simply the fraction of states with energies less than or equal to X(n, j). Note that 
X(n,j) and X(X(n, j), j) are, with respect to their first arguments, essentially inverses 
of each other. We again seek asymptotic expansions of these functions in the form 

X(X(n,j),j) = l(°\X(n,j))+lW(X(n,j))/j + 0(j- 2 ) (51) 
A(n, j) = \®(n/D) + X (l) (n/D)/j + 0(r 2 ) (52) 

where D = 2j + 1. Let {nj} be a sequence of labels such that 
limj^oo nj/ (2j + 1) = x G [0, 1]. We then define 

A (0) (x) = lim XOrijJ) and X {1) {x) = lim j \X{n v j) - X (0 \nJ{2j + 1))1 . (53) 

j— >oo j— >oo 

Similarly, if {nj} is such that linx^oo X(nj, j) = X define 

n.j 
2J + 1 

As to the existence and uniqueness of these limits we can say the following. First, it 
is clear that I^(A) equals T>^°\X) and that X^°'(x) is simply its inverse. A necessary 
condition for the existence of X^(x) is then that 

lim j [X(n, + m,j) - X(n 3 ,j)} = m/(2p^(X)) (55) 

must hold for any integer m. For j sufficiently large the gap A nj - = X(n + 1, j) — X(n,j) 
between successive eigenvalues must therefore behave as A„j w (2jp^(X(n, j)))^ 1 and 
be a slowly varying function of n in the sense that (A n+lj - — A n ,-) = 0(j~ 2 ). In section 



JW(A) = lim and X (1) (A) = lim j 

j-toc 2j + 1 j->oo 



I (0) (A(» 3 ,J)) 



(54) 
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I6.3l we consider a matrix sequence for which this condition is met only in certain regions 
of the spectrum. We will then also postulate a sufficient condition for the existence of 
\W(x) and XW(x). 



Assuming that the expansions in (1511) and (1521) are valid we proceed to investigate 
the relation between Z^(A) and X>W(A). First note that 



n i 



tr 



lim 

2j + 1 7 ^o+ 2j + 1 

1 + tanh(x/7))/2 and where 1/2 has been added to compensate for the 
fact that 7 (O) = 1/2. Combining this representation of X(X(rij,j),j) with ( 15TI) . ( 1521) 
and the results of section [5] leads to 



e 7 (A(n„j)-4) +1/2 



(56) 



with 6 7 (a; 



(A) = (A) and X (1) (A) = X? (1) (A) + ^ 



(57) 



for A G [A_,A+]. This agrees with the expectation that the eigenvalue distribution 
function is almost exactly the fraction of states with energies at or below A. The 
significance of the additional 1/4 term will become clear when we compare these results 
to numerical data. In cases where Z^(A) does not exist in the sense of fl54|) we may 
interpret I?W(A) + 1/4 as representing X^'(A) in an average sense. 



6.2. Expectation values and eigenstates 

In this section we work to lowest order in 1/j. Consider a sequence of tridiagonal 
Hermitian matrices {Qj} defined by the functions Aq(x) and Bq(x) as in ( 1131) . 
Let \X(n,j)) denote the eigenstate of Hj with eigenvalue X(n,j) and expansion 
coefficients {c m = (j, — j + m\X(n, j))}^ =0 in the Jz-basis. The expectation value 
(A(n, j)\Qj\X(n, j)) can then be expressed as 



(X{n,j)\Qj\X(n, j)) = / dxA (x)ip nd (x) + 2 / dxB (x)<f) nd (x 



where 



X 



2i 



m=0 



X 



m 
2] 



and 



2i-i 

m=0 



X — 



2m + 1 

It - 



(58) 



(59) 



These functions encode information regarding the norm and relative signs of the 
eigenstate's expansion coefficients. Let {nj} be a sequence of state labels such that 
lim^oo X(rij,j) = X. We will assume that the weak limits if)\(x) and 4>\{x) of i/j n .j(x) 
and <p njt j(x) exist and are determined by A alone. The same will then hold for the 
expectation value, and we may define 



(Q)(A)= lim (\(nj,j)\Qj\\(nj,j)). 



(60) 
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0.0 



0.5 



Figure 2. (a) a(x) and /3(x) for the example considered in section [6731 The form 
of il>\( x ) at A = —8/5, —1/2, 1/2 is also shown, (b) Numerical expectation values for 
Qj = J z /(2j) + 1/2 with respect to the eigenstates \X(n,j)) of Hj. Here j = 30. The 
solid line is the result of (|63[) . The dashed lines are obtained from ((53"]) by restricting 
the integration domain to one of the intervals that define the support of ip\(x). 



The "wave functions" ip\{x) and <f)\{x) can be related to the asymptotic distribution of 
eigenvalues as follows. Consider the matrix sequence {Hj + eQj} and denote the n'th 
eigenvalue of Hj + eQj by X e (n,j). The Hellmann-Feynman theorem now implies that 

dX e (n,j) 



(X(n,j)\Qj\X(n,j)) 



de 



= /dxioW^M + /dxBoW^. (61) 
e=0 J oA (x) J oB (x) 

Working to lowest order in l/j we replace X(n,j) by X^°'(n/D) and note that from (}5Tj) . 
(152]) and §7]) we have n/D = X^A^ra/D)) = Z>(°)(A( >(n/.D)). Taking functional 
derivatives to A (x) on both sides of the latter equation then leads to 

5\V*(n/D) -1 5V^(X) 



(62) 

\=\(0)(n/D) 



SA (x) p(°)(A) <L4 (a; 
with a similar result for -Bq(x). Performing the functional derivatives produce 

(<5>(A) = y dxi (x)^ A (x) + 2 y dxB (x)Mx) (63) 

where 

1 A — A (x) 

M x ) = ^ff^V [a{xmx)] {X) and (f> x (x) = 2Eo(x)p (o) (A) ^K^)^)]( A )- ( 64 ) 

with f2[ Qi/ a](A) as defined in (|48|) . Note that if>x(x) and 0a (%) are supported on the set 
S(X) = {x : a(x) < A < /%)}. 

Example: Alternating States 

A general sufficient condition for X^(A) and (Q)(X) to be well-defined in terms of the 
limits in (|54p and flBDl) is still lacking. In this section we investigate this question in the 
context of a specific example. The insight gained from this special case will suggest a 
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simple condition under which the results of the previous two sections may be applicable. 

Consider the sequence defined by A (x) = — 10x 2 + llx — 5/2, B (x) = x(l — x), 
Ai(x) = and Bx(x) = 0. It was found that the functions ij)\(x) and (f>\{x) are supported 
on S(X) and encode information regarding the asymptotic structure of the eigenstates. 
The form of S(X) therefore reveals how the states are localized in the J z basis. Figure 
[2(a) shows a(x) and (3(x) as well as the shape of <fi\(x) at three different energies. (See 
also figure [3] (b) for a comparison with numerical results for the Lipkin model.) The 
spectrum can be divided into three regions based on the number of solutions to a(x) = X 
and (3(x) = A: (i) A < -3/2, (ii) -3/2 < A < 1/32 and (iii) 1/32 < A < 1 where 1/32 is 
the maximum value of oc(x). In regions (i) and (iii) we see that S(X) is a single interval 
while in (ii) it is the disjoint union of two intervals. Since this sequence is closed these 
intervals are always terminated by solutions to a(x) = A and {3(x) = A. 

An important fact which is not reflected by the asymptotic results is that at finite 
j the eigenstates in region (ii) are localised alternatingly in the two intervals which form 
<S(A), but never in both simultaneously. This detail is lost in the large j-limit which 
effectively averages the properties of neighbouring eigenstates. This is also seen in the 
numeric results for the expectation value of Qj = J z / (2j) + 1/2 shown in figure [2] (b). In 
regions (i) and (iii) the expectation value varies slowly as a function of \(n,j) but alter- 
nate between two branches in region (ii). The prediction of fl63l) and fl64|) is indicated by 
a solid line and, while agreeing with the numerical results in regions (i) and (iii), clearly 
produces an average value for (X(n, j)\Qj\X(n, j)) in region (ii). The correct expectation 
value along a single branch may be found by restricting the integral in (|63l) to one of the 
two intervals that form S(X). To preserve normalisation the same restriction is required 
in the integral for p(°)(A). The results of this procedure appear as dashed lines in the 
figure. Similarly, we note that in region (ii) the gap A n j = A(n + 1, j) — A(n, j) also 
alternates between two branches and is therefore not a slowly varying function of n as 
required by f l55l . 

We are led to conclude that in region (ii) the quantities X^(A) and (Q)(X) are not 
well-defined in terms of the limits in (154"]) and ( l60l) since the results will depend on the 
choice of the index sequence {%}. In this region the analytic expressions for Z^(A) and 
(Q)(X) in ( |571) and (1631) must be interpreted as representing averages over neighbouring 
states in the spectrum. Based on these observations we postulate that ZW(A) and 
(Q)(X) are well-defined according to ( |54l) and ( |63l) for regions of the spectrum where 
iS(A) is a single interval. This phenomenon was also noted in [T4"l [15] in the context of 
the Lipkin model. 
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The integral expressions for D(°)(A) and 2?W(A) in flU} and (SSI) can easily be generalized 
to apply to closed sequences of Hermitian matrices with a fixed number M of smooth 
off-diagonal bands. The matrix elements, from top left to bottom right, of the m'th 
lower band is denoted 



(m) 
l n,j 



c {m) + i d 



(m) 



for 



where d. 



(0) 

71. J 



0. The matrix structure is 



Hi 



r (°) 


sP> 


a? 


4 ] 






(2) 


4" 


(0) 


a (M) 









(Af) 
a M+l 













n — m, 



2j 



(65) 



*2j 



;(Af) 
'M+l 



;(Af) 
l 2j 



s (l) 



(7 



l 2j 



(7 



2j 

(!) „(°) 
2j 



(66) 



2j J 

5 M - 



where, for compactness, we have dropped the j subscript and use ah, ' = ( a 
smoothness condition in ffT3l) now becomes 



(m)> 



S m ) 


= ct ] ( 


' 2n — m 


K 4j 




n( m ) 
— u 


/2ra — m 
V 4j 



and closure requires that Cq" 1 ^) = Dq ~'{x) 
Tq(x,9) and 0) by 

Af 

jr fc ( X; = cf \x) + 2 ^ fe m) (x) cos(m#) + D { ™\x) sin(m#) 



4- -r^ 



(m) 



4j 
2n — m 



4j 



k . The 

(67) 
(68) 



for x = 0, 1 and m > 0. Now define 



(69) 



m=l 



Repeating the arguments of section [3] leads us to conclude that H^°\x,9) = Tq{x, 
while 

I / „v (1 — x)x , 



-9|j* (x, 



-c^FoM). (70) 



16(1 - x)i P " " v ' 4 
Since the sequence is closed there are no edge corrections and the arguments of sections 
H] and [5] can be employed virtually unchanged to arrive at general expressions for the 
lowest and linear order contributions to the asymptotic eigenvalue distribution: 



Z><°>(A) 



dx 



2 Ti- 



de 

2^ 



eCA-JbM)) 



;>1 /-27T jfl 

2?W(A)= y dx ^ _5(A-J- (x,e)) 



Ji(x,0) + ^^3 a JbM) 



(71) 
(72) 
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Taking functional derivatives to Cg (x) and (x) now produces integral expressions 
for expectation values analogous to ( 1631) and ( 1641 ). 



8. Application I: Collective Spin models 

Consider a collection of N spin- 1/2 systems interacting via an infinitely long-ranged 
interaction in the presence of an external magnetic field. Since all the spin pairs interact 
identically the Hamiltonian can be expressed in terms of collective degrees freedom which 
are just the three components of the total spin 5*. The most general Hamiltonian of this 
type is [16] 

H = h ■ S + gj 2 x + g y S 2 y . (73) 

This model remains of great interest, particularly because it allows for the study of 
non-trivial quantum critical phenomenon in a simple setting. Here we will demonstrate 
how such models can be treated within our formalism in a simple and straightforward 
manner. For discussions of the physical phenomenon see [161 EH E] an d references 
therein. 



In the ferromagnetic phase where g x ,g y < the ground state belongs to the 
subspace on which S 2 = s(s — 1) with s = N/2. Focusing on this irrep of su(2) 
the matrix representation of H in the S^-basis is seen to have at most two off-diagonal 
bands with matrix elements exhibiting the smoothness property of ( 1681) with j = s. 
Since the latter is inherited directly from the su{2) generators it is clear that any 
Hamiltonian constructed in terms of these generators can be treated in exactly the same 
way. Furthermore, these sequences are always closed and the symbols can be calculated 
exactly, even at finite j, by using the representation of the generators as differential 
operators [12]. After rescaling g x>y — > g x , y /N and H — > H/j we find according to (169|) 

^(x, 0) = (2x - l)h z + (1 - x)x (g x + g y ) - 2a/ (1 - x)xh y sm(9) 

+ 2a/(1 - x)xh x cos(0) + (1 - x)x (g x - g y ) cos(26») (74) 
1 



/ x h x cos(9) — hy sin(9) 

g x + g y + (g x -g y )cos(29) + - x 



(75) 



a/(1 -x)x 

Inserting these expressions into ( J72l) yields explicit integral expressions for the eigenvalue 
distribution function and its finite size corrections. Taking derivatives to the various 
coupling constants then produce expressions for expectation values, although in this 
regard the caveat discussed in section 16.31 must be kept in mind. We will not embark 
on a study of the general Hamiltonian here, but instead highlight two prominent special 
cases. We return to the notation used for tridiagonal sequences prior to section [71 
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8.1. Lipkin Model 

As a special case of the well-known LGM model [T8] . itself a particular instance of (!73|) . 
we consider 

^ = 4 + f(^ + 5 2 ) (76) 

within the s = iV/2 sector. Note that H leaves the even H e = span{|s, — s + n) : n even} 
and odd "H a = span{|s, — s + n) : n odd} subspaces invariant, and that its matrix rep- 
resentation in each of these is tridiagonal. We define j = s/2 and j = {s — l)/2 in the 
even and odd sector respectively. After rescaling the matrix elements by 2j to render 
them extensive we find that, in both sectors, Aq(x) — 2x — 1, Bq(x) — A(l — x)x and 
Ai(x) = 0. The two sectors are therefore distinguished by Bi(x) which is Bi(x) = A/8 
and Bi(x) = A(3 — Ax + 4x 2 )/8 for the even and odd case respectively. Deriving the 
asymptotic distribution function X>(°)(A) and its finite-size correction V^(X) is now 
simply a matter of evaluating the integrals in (j4"2"]) and fj4Tl) . This can be done exactly 
with the final result given in terms of elliptic integrals. The results agree with those 
obtained in [TJJ [15] based on an analysis of the zeros of the Majorana representation of 
the eigenstates.We will not analyse these results further here, but only add the following 
comments. 

Referring back to f j48l) and ( |64l) we see that since (A — a(x))((3(x) — A) is a fourth 
degree polynomial in x the calculation of expectations values for observables defined 
by polynomial A (x) and B (x) using (1631) will involve evaluating elliptic integrals. 
These integrals are known to satisfy certain recurrence relations [19] which aid in 
this calculation. For example, consider the observable Qj = J z /(2j) + 1/2 for which 
Aq(x) = x and Bq(x) = 0. The moments of Q then satisfy 

(<T>(A) = 2a * _ 1} ^{2 + 1- 2m) (g m -)(A) (77) 

^ ' i=l 

where (a , a h a 2 , a 3 ) = (47 2 , -87 2 , 4 (7 2 - 1) , 4(A + 1), -(A + l) 2 ). This relation is ex- 
act in the thermodynamic limit and holds for all energies and values of the coupling 
constant. Once (Q)(X) and (Q 2 )(X) are known from ( 1631) all the higher moments of Q 
therefore follow recursively from (IT7j) . 

Finally we perform some comparisons with numeric results. LetjW(A) and£>( L )(A) 
denote the linear order (in 1/j) approximations to X(A) and T>(\). Keep in mind that 
from (J57]) it holds that X (i) (A) = V^ L \\) + l/(4j). In figure [3] (a) the three smooth 
curves, from top to bottom, are X^ L ^(A), T>^ L \\) and T)( L \\) — l/(4j). The dots indicate 
exact numerical values for j = 60 while the staircase curve S(X) indicates the fraction 
of eigenvalues at or below a certain energy. The three smooth curves clearly represent 
slightly different large-j approximations to S(X). In particular, we see why X(A), rather 
than T>(\), is the true reflection of the number of eigenvalues less than or equal to a 
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Figure 3. (a) A comparison of analytic and numerical results for j = 60 and 
7 = 2.5. See text for details, (b) Plots of the functions ip\(x) and 4>\(x) for the 
Lipkin model with 7 = 2.5 and A s» —0.67. Dots correspond to numerical data for 
j = 1500. In the main plot the dots appear at (m/(2j),2jc^ n ) and in the inset at 
((2m — l)/(4j), 2jc m c m +i). The dots therefore represent the magnitudes of the 6- 
functions in ([59]) . but scaled by 2j to aid comparison with ip\(x) and 4>\(x). 



particular X(n,j). Figure [3] (b) compares numerical results for ip n j{x) and (f> n> j(x) at 
j = 1500 with ij)\{x) and 4>\{x). 

8.2. Uniaxial Model 

Another prominent special case of f!73|) is the uniaxial model for which the Hamiltonian, 
in the symmetric phase [17] and after appropriate rotations, reads 

H = S X - ^S 2 Z . (78) 

Setting j — s — N/2 and rescaling the H by j to render it intensive then pro- 
duces the tridiagonal matrix sequence {Hj}- The symbols can be calculated easily 
by representing the su{2) generators as differential operators [12] , and we find that 
A (x) = -27(1 - 2x) 2 , B (x) = y/(x-l)x, Ax{x) = and B x {x) = 1/^64(1 — x)x. A 
myriad of information concerning the spectrum and expectation values, in both phases 
and for all energies, now follows from a straightforward application of the results derived 
in sections [5] and El Here we only present new analytic results for the asymptotic density 
of states. 

Using (144 j) it is found that in the thermodynamic limit the eigenvalues of Hj range 

from 

A (7) = I _1 7 K 1/4 (79) 

Al) { -27-1/(87) 7 >V4 U9J 

up to A + = 1. We first consider the energy range A G [—1,1]. In terms of 
A = (1 + 87A + 167 2 ) 1/2 and B = (A - 1 - 4r/\)/(2A) the density of states in the 
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thermodynamic limit is p(°)(A) = K(B)/(\/Att) while its finite size correction reads 
(1) = S(X-l) *(A + 1) d [2^W-(l + J 4 + 8A 7 )A- W l 

Here K(x) = F (fix) and E(x) = E |x) are complete elliptic integrals of the first 



and second kind 

Now consider A G [A_(7), -1] with 7 > 1/4. We set C± = -1 - 4A7 ± 47VA 2 - 1 
and find that p(°>(A) = 2K(C_/C + )/(tt v CV) and 

~C+E(C-/C + ) + (C_ + \)K{C-/C+ 



P ( 1 )(A) = 2 ^ (A - A - ) + A 

7I67 2 - 1 dA |_ 47r 7v ^ 
9. Application II: Orthogonal Polynomials 



.(81) 



It is well known that any sequence of orthogonal polynomials {p m (x)}™ =0 satisfies a 
three-term recurrence relation 1 20 1 of the form 



xp m {x) = b m+1 p m+ i(x) + a m p m (x) + b m p m -i(x) (82) 

where b m > 0. The symmetric, tridiagonal Jacobi matrbjj] J m is then defined in 
terms of the recurrence coefficients as [J m ] n ,n = <x n -i an d [J m }n+i,n = [J m }n,n+i = b n 
for n — 1, . . . ,m. Its eigenvalues are precisely the m zeroes of p m (x). This simple 
correspondence enables the calculation of the asymptotic zero distributions plus 
corrections for a range of orthogonal polynomial sequences. Here we present only two 
well-known cases. Examples of other polynomial sequences to which this formalism can 
be applied may be found in [SJ [20] . 

9.1. Laguerre Polynomials 

The recurrence coefficients for the generalized Laguerre polynomial L^x) are 
a n = 2n + a + 1 and b n = yn(n + a) where a > —1 is a real parameter. We define 
the matrix sequence {Hj} by Hj = J 2 j+i/j. The n'th eigenvalue X(n,j) of Hj and 
root v(n,j) of Lf,- +1 (rr) are therefore related by X(n,j) = b>(n,j)/j. We also allow a 
to depend on j as a = aoj + «i- The eigenvalues of Hj and the zeros of L^+i*" 
therefore share the same asymptotic distribution. The matrix elements of Hj are 

a nJ = (2n + a j + a 1 + l)/j and b nJ = y / n(n + a j + ai)/j (83) 

from which it follows by ([TBI that 

A (x) = 4x + a 5o(.t) = a/22 (2x + a ) (84) 

= 1 + ax Sx(x) = —7=== f (85) 

4a/2x (2x + a ) 

Note that this sequence is not closed on the right since B (l) > 0. 



I Not to be confused with the su(2) generators. 
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From this we can derive analytic expressions for Z>(°)(A) and PW(A) which appear 
in the appendix. The corresponding corrections to the density of zeros are 

P (°)(A) = ^HMZZ3 (86) 

(1) = (A - a ) (A - a + 2 Ql - 2) - 2a _ 1 _ 
8ttAv/(A + -A)(A-A_) 8 

- h(X - A_) (1 + 2a 1( 5 QO ,o) (87) 



which are supported on A G [A_,A + ] with A± = 4 + a ± 2a/4 + 2«o- Here 5 a) & is 
the Kronecker- delta. These results agree with those found in [2TJ for a closely related 
sequence. These functions are depicted in figure H] (a) for a specific choice of parameters. 
As a concrete measure of the accuracy of the asymptotic expansion we consider the 
problem of finding the n'th root of L^^ ai \x) for a given large j. This involves solving 
for X(n,j) numerically from 

I^(X(n,j)) + tr«(A(n,j)) 



j v v 2 j + l 

where 1^(X) is related to D^A) by ([57}. The errors in the resulting approximations 
to X{n,j) appear in figure H] (b) as a function of n for j = 250, «o = 1 and ct\ = 5. 
From curve (i) we see very good agreement with the exact results with errors being at 
worst ~ 0.2% and on average only « 0.005%. Using the lowest order approximation (i.e. 
solving 1^ (X(n, j)) = n/(2j + 1)) produces errors of about two orders of magnitude 
greater (curve (ii)); in line with the fact that j = 250. An insighful alternative measure of 
the accuracy is to compare the error in the approximation of X(n,j) with the difference 
between X(n,j) and X(n + 1, j). This gives a indication of the "resolving power" of 
this method, i.e. whether the approximation of a single root is sufficiently accurate to 
reliably distinguish it from its neigbours. When working up to lowest order this is not 
possible, as indicated by curve (iv); the error is of the same magnitude as the distance 
between succesive roots. However, to linear order (curve (iii)) we find that the errors 
are at most ~ 2.7% and on average only w 0.25% of this distance. Since the expressions 
for Z( 0,1 )(A) involve only elementary functions this procedure amounts to a simple and 
efficient numerical algorithm for accurately approximating arbitrary roots of high degree 
Laguerre polynomials. 



9.2. Jacobi Polynomials 



The recurrence coefficients for the Jacobi polynomials Jn^\x) are 

P 2 -a 2 

P) = (2n + a + m + 2n + a + P) (89) 

hi R\ - 9 / n(n + a)(n + f3)(n + a + f3) 

o n{ a, P) ^ ^ a + f3 + 2n _ l)[2n + (y + ma + f3 + 2n + 1) W 

where a > — 1 and > — 1 are real parameters. We again allow these parameters to 
vary with j as a = a j + a\ and (3 = (3 j + f5\. The precise polynomial sequence is 
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Figure 4. Results for the Laguerre polynomials with — 1 and a.\ — 5. (a) The 
lowest and linear order contributions to the asymptotic density of zeros, (b) Errors 
in the results of approximating X(n,j) by solving (|88[) for a given n when j = 250. 
Interpolating curves rather than individual points are shown. Solid curves (i) and (ii) 
show the percentage error when using the linear and lowest order approximations to 
2j(A) respectively. Curve (iii) shows the error at linear order expressed as a percentage 
of A(n + 1, j) — A(n, j). Curve (iv) is as (iii) but for the lowest order approximation. 



therefore J^^^'^^^lx) while the corresponding matrix sequence {Hj = J2j+i} has 
elements 

a n ,j = a n (a j + cm, /3 j + ft) and b nJ = b n (a j + a x , (3 j + ft). (91) 

Expressions for {Ao,i,5o,i} follow easily from ( fl3i) and ( 1911 . Performing the integrals 
yields 



ro) m _ («o + ft + 4)y/(A + -A)(A-A_) , 
P [A) ~ 47rfl - A 2 > 

p«(A) 



(i)/v, a + ft + 4 



4ttv/(A + -A)(A-A_) 



a a ft/3 (a + ^-2) 



1 - A 1 + A 2(a + ft + 4) 



- ^5(A - A+) [1 + 2ft5 A)i0 ] - l5(X - A_) [1 + 2a 1( 5 Qo ,o] (93) 

o o 

where a = «o — 2«i and /3 = ft — 2ft. These densities are supported on A G [A_, A + ] 
with \± = Aq(1) ± 2Bq(1). These expressions differ from those obtained [22j for the 
same polynomial sequence. We have performed numerous numerical tests to verify our 
results and have not found any discrepancies. As in the case of the Laguerre poly- 
nomials analytic expressions for Z^'^A) can be also be found. These allow for the 
implementation of the same numerical algorithm for finding specific roots of high degree 
Jacobi polynomials. Comparison with exact results show the same levels of accuracy as 
observed in the Laguerre case. 
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10. Conclusion 

We have presented an approach to studying the asymptotic properties of a class of 
matrix sequences which arise in a number of diverse settings. Analytic expressions for 
the asymptotic eigenvalue distribution, eigenstates and expectation values were derived. 
Underpinning this approach is idea of reformulating matrix (or operator) problems in 
terms of their associated symbols for which the multiplication rule is given by the star 
product. A systematic expansion of the star product in orders of the inverse matrix 
dimension (1/D) then allowed for the derivation of a general trace formula which cap- 
tures both the limiting D — > oo behaviour as well as linear order corrections. Once a 
particular matrix sequence has been specified through the functions {A 0> i(x), B 0> i(x)} 
applying these results simply amounts to evaluating certain integral expressions. We 
considered applications to quantum spin systems and orthogonal polynomials and found 
that in both cases new and existing results can be obtained in a straightforward manner. 
In fact, in these two fields alone a range of applications are still to be investigated. A 
thorough investigation into applications in other fields should also prove fruitful. 

Two aspects of our approach can benefit from further generalization. First, one 
would like to relax the closure restriction on the band diagonal case considered in section 
[3 This will require a more general analysis of the edge corrections. Existing results 
regarding Toeplitz matrices may be of use in this regard. Secondly, one may consider 
the inclusion of higher order corrections going beyond the linear case. The existence of 
the star product expansion to all order, as reflected in expression (JS}, suggests that this 
is certainly a possibility. 
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Appendix A. Edge corrections 

Here we derive expressions for the edge terms 7^ and 7^ appearing in the linear order 
correction to the trace. Focussing on 72 we see from ( 1251) and (T57j) that 

= lim lim f dx f ——j 

Note that the ^-integral acts as a projection which eliminates off-diagonal contributions 
in both (z\f(Hj)\z) and f(H^). Only the diagonal entries of f(Hj) are relevant, 
and similarly only the first term in the expansion f(H^) = ^ n=0 9n{ x ) cos(n9) 
contributes. With this in mind we define Dj as the diagonal part of f{Hj) and 



(A.l) 
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D(°\x) = (27T)- 1 d6 f(H^) in terms of which lf ] reads 

lf ] = lim lim / dx j \(z\Dj\z) - D {0) (x)] (A.2) 

e->0+ j^oo J L 



with z = a/x/(1 — x). 



We first demonstrate the calculation of 7^ for the case where Aq(0) = and 
f(x) = x 4 . Let do, ... , d 2 j denote the diagonal matrix elements of Dj = diag(ff 4 ). Up 
to lowest order in 1/j it holds that 

d =a A + Qa 2 bl + 2b\ + 0(1/ j) (A.3) 
d x = a\ + Yla\b\ + 56 4 + C(l/j) (A.4) 
d n = a 4 n + 12a 2 n b 2 n+ i + 66 4 +1 + 0(1/ j) 2 < n < 2j - 2 (A.5) 

where the smoothness condition in ( Tl3|) has been used to simplify the expressions by 
replacing, for example, a® + a,\ by 2ao + 0(1/ j) and so on. Due to the tridiagonal 
structure of Hj the numeric coefficient of the b^ +l term in d n has a simple combinatoric 
interpretation: it is the number of four step random walks on the lattice {1,2,3, . . .} 
which start and end at site Since only the walks associated with c? an d c?i 

are affected by the lattice boundary the numeric coefficients appearing in do and d\ 
do not conform to the general pattern observed in both d n for 2 < n < 2j — 2 and 
in D^(x) = A (x) 4 + 12Ao(x) 2 B (x) 2 + 6B Q (x) A . (The coefficient of the a 2 b\ term is 
similarly affected, although this will be seen not to influence 7^ .) It is this deviation 
of the numeric coefficients in do and d\ from those in D^(x) that gives rise to non-zero 
edge corrections. To make this explicit, let Dj denote the diagonal matrix sequence in 
which these deviations in do and d\ have been "corrected" to lowest order by defining 
d' n = d n for 2 < n < 2j — 2 and d' Q = do + 46 4 and d[ = do + b\. Note that no corrections 
are necessary for terms containing a n since Ao(0) = implies that a n = 0(1/ j) if 
n = O(j ). The matrix elements of D'- now represent, to lowest order in 1/j, a sampling 
of the function D^(x) at the points n/(2j) with n = 0, . . . , 2j. It is clear that both {D'j} 
and {Dj} converge to D^°'(x) on x G (0, 1) but it also holds that j[(z\D'j\z) - D^(x)} 
remains bounded as j — > oo. Replacing Dj by in ( 1A.2I) therefore allows the j — > oo 
limit to be taken into the integral after which taking e — ?■ + produces a zero result. It 
follows that 

TjV = lim lim / dxj{z\Dj - D'Az). (A.6) 

Furthermore, since only the first two diagonal matrix elements of Dj — Dj are nonzero 
it holds that j(z\Dj — D'j\z) will become completely localised around x = as j — > oo 
and so lim^oo f* dx j(z\Dj — D'j\z) is independent of e G (0,1). Combining these 
observations with the expression for the trace in (jl]) leads to 

7 L (1) = lim lim / dxj(z\Dj -D'Az) = lim - V(d nj -d' nj ) (A.7) 

e-)-0+ i^oo J Q J j->oo 2 ^—^ 

u n=0 
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where 77 G (0, 1) is an arbitrary constant. The right hand side is readily evaluated to 
produce the left edge correction 

(A.8) 

x m with m even. Arguing as 



7^ = lim (-46? - bt)/2 = -5(B (0)) 4 /2. 



j-s-oo 



This result can be generalised to apply to any f(x) 
before the matrix elements of D'a now become 

< n < m/2 



d' 



d n + F m ^ +1 < n < m/2 (A.9) 
d' n = d n m/2 < n (A.10) 
where F m ^ n is the number of m-step random walks on Z which start and end at n + 1 and 
visit site zero at least once. A standard argument gives F mn = ( n+ ™ / 2 ) anc ^ ex P ress i° n 
(1 A. 71) produces the left edge correction for f(x 

(Bo(0)) m ^ p _ (B (0)) 



x as 



r; 



(i) 



m 
m/2 



(All) 



2 / j '">" 1 

n=0 

whenever m is even. If m is odd any correction terms added to d n to match the pattern 
in D(°'(x) will contain at least one a, factor and are therefore zero in the j — > 00 
limit when A (0) = 0. 7^ therefore vanishes whenever /(x) is a odd power. For a 
general analytic f(x) expressed as a power series the edge correction will be a linear 
combination of the expressions in flA.llj) . It would be preferable to have an expression 
for directly in terms of f(x) itself. For this purpose we note that 



dx (2x) 



6{x-l) 5{x + l) 



,-kVT— X" ^ ^ 

and comparison with ( IA.11I) then reveals that, for a general analytic f(x), 

" +1 dx /(2Bq(0)x) 

-.1 ' 471^1 - x 2 

The only remaining restriction on this result is that Ao(0) must be zero, but this can be 
circumvented as follows. If a given {Hj} does not satisfy A (0) = then the sequence 
{Hj — v4o(0)} surely does. Instead of applying (1A.13I) to the sequence {Hj} and function 
f(x) we simply consider {Hj — A (0)} and g(x) = f(x + A o (0)) instead. Recalling that 
a(x) = A (x) — 2B (x) and (5{x) = A (x) + 2B (x) the final expression for can be 
written as 







( m ) - T 
Km/ 2) L 



m odd 
m even 



(A.12) 



-/(-2So(0))/8-/(2B (0))/8 



(A.13) 



(i) 



r, 



where 



/(a(O))/8-/G0(O))/8 + 



/3(0) 



dx 



a(0) 47Ty/(P(0)-x)(x-a(0)) 



x) 



(47T) 





arcsm 



The corresponding expression for follows by replacing [a(0),/5(0)] with [a(l),/3(l)] 



Q+/3-2A 



A G [a, (3} 
otherwise 



(A.14) 
(A.15) 

(A.16) 



in flA~T5|) . Finally we note that if B (0) = (B (l) = 0) then 7£ ' (7^ ') is zero. In 
particular, if the sequence is closed there are no edge contributions to ■ 
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Appendix B. and for the Laguerre polynomials 

Defining A = (47r)- 1 y / 8A - (A - a f, 

B = i arccos f^±£] and C = J_ arccos ( /*° + A I (B.l) 



7T 



V2v/2ao + V 2tt I 2^^ + 2) A 



we find that p(°)(A) = A + B-a C and £>W(A) = -A/2 + B/4 + (a -2a 1 )C/2-l/8. 
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